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Simulated Energetic Particle Transport in the Interplanetary 
Space: The Palmer Consensus Revisited 

R. C. Tautz^ and A. Shalchi,^ 

Abstract. Reproducing measurements of the scattering mean free paths for energetic 
particles propagating through the solar system has been a major problem in space physics. 
The pioneering work of Bieber et al. [Astrophys. J. 420, 294 (1994)] provided a theo- 
retical explanation of such observations, which, however, was based on assumptions such 
as the questionable hypothesis that quasi-linear theory is correct for parallel diffusion. 
By employing a hybrid plasma- wave/magnetostatic turbulence model, a test-particle code 
is used to investigate the scattering of energetic particles. The results show excellent agree- 
ment with solar wind observations. 



1. Introduction 

The transport of energetic particles in a turbulent 
medium has been considered a key problem in various sub- 
fields of astrophysics. However, only in the solar system can 
one obtain detailed information about the transport parame- 
ters from measurements. Therefore, the solar wind provides 
the perfect laboratory to test our understanding of parti- 
cle propagation through a tenuous plasma. In his pioneer- 
ing paper, Palmer [1982] presented a detailed comparison 
of analytical particle diffusion coefficients along the mean 
magnetic field with real data. He found that the under- 
standing of particle transport theory was quite incomplete 
at this time. The theoretical description for this comparison 
was based on quasi-linear theory combined with a simplified 
turbulence model (magnetostatic slab without turbulence 
dissipation range). 

A few years later, Bteber et al. [1994] improved the the- 
oretical description of particle transport by incorporating 
three important effects: (i) the one-dimensional slab model 
was replaced by a quasi three-dimensional model known as 
the slab/2D composite model. In the two-component model, 
the turbulence is described as a superposition of slab modes 
and so-called two-dimensional (2D) modes; (ii) dynamical 
turbulence efi'ects were included to account for the fact that 
magnetic field vectors decorrelate with respect to the initial 
conditions as time passes; (iii) a strongly decreasing part 
of the turbulence spectrum was included, which is known 
as the dissipation range (see [Schlickeiser et al, 2010; Chen 
et al, 2010] and references therein). Therefore, a complete 
spectrum was used that contained energy, inertial, as well as 
dissipation ranges. As demonstrated by Bieber et al. [1994], 
solar wind observations of energetic particles were success- 
fully reproduced. 

Since then our understanding of plasma-particle interac- 
tions has significantly improved (see, e. g., [Shalchi, 2009] for 
a review) . It was shown that quasi- linear theory itself, which 
was the basis for [Bieber et al., 1994], is highly questionable 
[Tautz et al., 2006] mainly due to: (i) the well-know 90° 
scattering problem describes the inability of the quasi-linear 



approach to describe correctly the pitch-angle scattering at 
pitch-angles close to 90°, which problem has been known 
since the early seventies of the 20"^ century [Volk, 1973] 
and has been revisited more recently [Tautz et al., 2008; 
Shalchi et al., 2009]; (ii) 90° scattering has been explored in 
numerical work [Qm and Shalchi, 2009] showing that quasi- 
linear theory is not reliable in general case; (iii) the so-called 
the geometry problem is ascribed to the observation [Shalchi 
et al., 2004; Shalchi, 2007] that quasi- linear theory overesti- 
mates the parallel mean free path for the standard slab/2D 
model used in [Bieber et al., 1994]. All three problems are 
related to the applicability of quasi-linear theory to model 
the parallel transport of energetic particles in the solar sys- 
tem. 

In this paper, in contrast, we revisit the Palmer consensus 
by using advanced numerical simulations. In test-particle 
codes, particle diffusion coefficients can be computed with- 
out employing a transport theory such as quasi-linear the- 
ory. Specifically, the so-called Padian code [Tautz, 2010a] 
will be used to compute the parallel diffusion coefficient for 
a hybrid turbulence model consisting of both plasma waves 
and magnetostatic turbulence and for the dissipative turbu- 
lence spectrum used in [Bieber et al., 1994]. In Sec. 2 we 
briefly discuss the turbulence model used in our simulations 
and in Sec. 3 we describe the test-particle code. In Sec. 4 
we present the new numerical results as well as a comparison 
with observations. In Sec. 5 we conclude and summarize. 



2. The Turbulence Model 

The magnetic field scenario adopted in the present work 
can be approximated by a superposition of a homogeneous 
mean magnetic field and a stochastic component, i.e., by 
B{x) — BoSz + SB{x). The first contribution corresponds 
to the solar magnetic field on a sufficiently small scale. A set 
of Cartesian coordinates is chosen so that the mean magnetic 
field is aligned with the z direction. The second contribu- 
tion is more difficult to model. The stochastic contribution 
is generated by the solar wind turbulence and is described 
in the following paragraph. 
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2.1. The Magnetic Correlation Tensor 

The fundamental quantity in the theory of magnetic tur- 
bulence is the time-dependent correlation tensor in wave- 
vector space. The tensor is defined as 

Pe,^{k,t) = {5Be{k,t)5B*^{k,0)}, (1) 

where we used the ensemble average (...) and the wavevec- 
tor k. A common assumption in turbulence and diffusion 
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theory is that all tensor components have the same temporal 
behavior, i.e., Pern{k,t) = P im{k)r{k,t), where the mag- 
netostatic tensor has [Matthaeus and Smith, 1981; Schlick- 
eiser, 2002], for axisymmetric turbulence, the form 



P<m(fc) =G'(fe||,fcj 



(2) 



Furthermore, we introduced the so-called dynamical correla- 
tion function r(k, t), which is discussed below. The function 
G(fc||,fcx) describes the scale dependence of the magnetic 
fluctuation and is discussed in what follows. 

2.2. Two-component Turbulence 

In the solar system, the standard model for the mag- 
netic field fluctuations is the so-called slab/2D composite 
model. In this model, which is sometimes called the two- 
component model, it is assumed that the spectrum can bo 
approximated by a superposition of slab modes defined as 
SB{x) — 5B{z) as well as two-dimensional (2D) modes with 
5B{x) = SB{x,y). Physically, the former modes represent 
parallel Alfvenic fluctuations and the latter modes ideal- 
ize the quasi-2D structures including flibx tubes that can 
develop from the interactions of such waves [Belcher and 
Davis, 1971; Shebalin et al., 1983]. 

The two-component model can be confirmed by using ex- 
tensive analyses of solar wind data [Matthaeus et al, 1990: 
Narita et al., 2010]. Such observations indicate that ap- 
proximately 80% of the turbulent inertial range energy bud- 
get is in the two-dimensional fluctuations [Bieber et al., 
1996]. The rest of the magnetic energy is contained in 
the slab modes. More recent data analysis projects con- 
firmed this structure [Horbury et al., 2008; Turner et al., 
2011, 2012]. Furthermore, there is a strong support of the 
aforementioned turbulence model by numerical simulations 
[Matthaeus et al., 1996; Shaikh and Zank, 2007; Dmitruk 
and Matthaeus, 2009] and analytical treatments of turbu- 
lence [Zank and Matthaeus, 1993]. Additionally, the slab/2D 
model is very tractable in analytical studies of magnetic field 
lines and energetic particles [Matthaeus et al., 1995; Bieber 
et al, 1994; Shalchi et al, 2010]. 

In the model described here, the turbulence power spec- 
trum has the form 



G(A;||,A;x)=/^''(fe||) 



5{k^ 
k± 



k± 



(3) 



where the onc-dimcnsional spectra are g^^^^ for slab modes 
and for two-dimensional modes. Furthermore, parallel 
and perpendicular wavenumber components are denoted by 
and k±, respectively, and 5{x) is the Dirac delta distri- 
bution. 

2.3. The Turbulence Spectra 

For the two turbulence spectra the form proposed by 
Bieber et al [1994] is used, namely for the slab spectrum 



27r 

[l+ (fc[|4lab)']" 



s/2 



for fc|| ^ fcslab 

[l + (A;slab4lab)^] (~f^) ' ^^"^ ^ ^"''^'^ 



From the condition that the wavenumber-integrated corre- 
lation tensor be equal to the total magnetic field energy, the 
normalization function is found to be 



C{s) 



r{s/2) 



Table 1. Suniinary of the various parameters aRRiiriied for 
the test-particle sijimialions. 



Parameter description 


Symbol 


Value 


Inertial range spectral index 


s 


o/o 


Damping range spectral index 


P 


Q 
O 


Alfven speed 


VA 


33.5 km/s 


Slab bend-over scale 


^slab 


0.03 AU 


2D bend-over scale 


i2D 


0.03 AU 


Slab dissipation wavenumber 


^slab 


3 X 10^ (AU)-i 


2D dissipation wavenumber 


fe2D 


3 X 10« (AU)-i 


Mean magnetic field strength 


So 


4.12 nT 


Slab turbulence strength 




0.2 (5B2 


2D turbulence strength 


SB'2D 


0.8 <5B2 



with the gamma function r{z). In the spectrum used here, 
we introduce the slab bendover scale fsiab, the slab dissipa- 
tion wavenumber fcgiab, the inertial range spectral index s, 
and the dissipation range spectral index p. Furthermore we 
denote the total magnetic field strength of the slab modes as 
''Bsiab- For the two-dimensional spectrum, g^^ , we employ 
precisely the same form, namely 



2D/, N 

g (k±) 



«2D0i>2E 

s/2 



for k± ^ k2u 



j ^1 + (A;x£2d) 

I [1 + (fc2D^2D)'] (l^y , for fcx > feD 



(5) 



The parameters have the same physical meaning as those 
used in the slab spectrum, Eq. (4). 

Of course, there are alternative models (e. g., [Sridhar and 
Goldreich, 1994; Goldreich and Sridhar, 1995]) that take 
into account the different physical origins of one- and two- 
dimensional fiuctuations (cf. [Matthaeus et al, 2007]). Here, 
however, identical spectra for both directions are assumed 
for reasons of (i) simplicity and (ii) to enable future com- 
parisons with non-linear calculations that are only tractable 
analytically for as simple as possible spectral distributions. 

2.4. Dynamical Turbulence 

The dynamical correlation function [Schlickeiser, 2002; 
Shalchi et al, 2006] can be expressed as 



r{k, t) = cos (wt) e 



(6) 



20F r((s-l)/2) 



where the trigonometric and exponential factors describe 
oscillations related to wave-propagation effects and damp- 
ing effects, respectively. Previously different models have 
been discussed to approximate the parameter a. Exam- 
ples are different plasma wave damping models [Achatz 
et al, 1993a, b], the damping model of dynamical turbulence 
[Bieber et al, 1994], and the non-linear anisotropic model of 
dynamical turbulence [Shalchi et al, 2006]. 

Here, we choose a = as test-particle codes are still 
unable to take into account damping effects. In contrast, 
the recently developed Padian code [Tautz, 2010a] can be 
applied for undamped plasma wave turbulence (see, e.g., 
[Tautz, 2010b]). Therefore, we use the simple model of 
undamped parallel propagating shear Alfven waves corre- 
sponding to u = jvAk\\ and a = 0, where the Alfven 
speed i;a = Bo/y/4-Kp represents the speed of the propa- 
gating plasma wave. The parameter j is used to track the 
wave direction (j = -|-1 for forward- moving and j = —1 
for backward-moving Alfven waves with respect to the am- 
bient magnetic field). Several studies have addressed the 
direction of propagation of Alfvenic turbulence [Bavassano, 
2003). In general it is expected that, closer to the Sun, most 
waves should propagate forward whereas at radial distances 
around 1 AU both wave intensities should be equal. In the 
current paper, therefore, we adopt the latter scenario. 
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Figure 1. (Color online) The parallel mean free path as a function of the particle rigidity for electrons 
(red diamonds) and for protons (blue squares) from the numerical simulations using the Padian code. 
The dissipation range spectral index is p = 3 and the dissipation wavenumber is fctj = 3 x 10^ (AU)~^. 
The comparison with the observational data (Palmer consensus range [Palmer, 1982]) shows excellent 
agreement. Data points obtained from Fig. 1 of Bieber et al. [1994], reproduced by permission of the 
American Astronomical Society (AAS). 



3. The Padian Code 

In the following we use a numerical Monte-Carlo code to 
compute the parallel diffusion coefficient of energetic parti- 
cles for the turbulence model described above. A general 
description of the code and the underlying numerical tech- 
niques can be found elsewhere [Tautz, 2010a]. Specifically, 
the generation of turbulent magnetic fields corresponding to 
the composite model proceeds as follows: The slab and the 
two-dimensional components are generated separately using 
[Tautz, 2010b] 

5B{r, t) = ?ltY2 e'±A{kn) exp{ i [k„z' + /3„ - uj{K)t] } , 

(7) 

where the wavenumbers k„ are distributed logarithmically 
in the interval fcmin ^ fcn ^ fcmax and where /3 is a random 
phase angle. The important part is the time-dependence, 
which is taken as (jj{k„) = jvA{kn)\\ for the slab modes (with 
j = 1 for n even and j ~ —1 for n odd) and io{kn) ~ for 
the two-dimensional modes. For the amplit ude and the po- 
larization vector, one has A{kn) oc \J g(kn) and • = 0, 
respectively, with the primed coordinates determined by a 
rotation matrix with random angles so that fe || for slab 
modes and fc _L for 2D modes. 

From the integration of the equation of motion, the par- 
allel diffusion coefficient, k\\ , is then calculated by averaging 
over an ensemble of particles and by determining the mean 
square displacement in the direction parallel to the back- 
ground magnetic field as Ky = (i;/3)A|| = {(Z\2:)^)/(2t). 

For the simulation runs, the various parameters are sum- 
marized in Table 1. For the minimum and maximum 
wavenumbers included in the turbulence generator, the fol- 
lowing considerations apply: (i) the resonance condition 
states that there has to be a wavenumber k so that Ri,k ^ 1, 



where i?L denotes the particle's Larmor radius so that scat- 
tering predominantly occurs when a particle can interact 
with a wave mode over a full gyration cycle; (ii) the scal- 
ing condition requires that Ri^Q^cit < imax, where Lma.^^ oc 
1/fcmin is the maximum size of the system, which is given by 
the lowest wavenumber (for which one has fcmin = 27r/Aniax, 
thereby proving the argument). In practice, the second con- 
dition determines the minimum wavenumber while the first 
condition determines the maximum wavenumber. 

Additionally, the relative strength of the turbulent mag- 
netic field as compared to the homogenenous field is as- 
sumed to be SB /Bo — 0.5 as suggested by Ruffolo et al. 
[2012]. Similar values were applied by He and Wan [2012] 
for the interplanetary medium near Earth. Analytical esti- 
mates based on a WKB approach [Zank et al., 1996] indicate 
that the relative turbulence strength is somewhat less than 
unity. 

Additional simulation parameters are chosen as fcmin^siab = 
10~* and fcmax^siab ~ 10 . Furthermore, the maximum sim- 
ulation time can be expressed as vt/lsish ~ 10^ and the sum 
in Eq. (7) extends over Nm ~ 2048 wave modes. 

An alternative method to simulate the transport of 
charged particles (so-called passive tracers) by following the 
evolution of the magnetic field in a magnetohydrodynamic 
description [Muller and Busse, 2007; Zank et at, 1996]. The 
comparison with either the Kolmogorov-based turbulence 
introduced above and with transport theories [Fumtko and 
Tooru, 2004] generally shows agreement; however, the dis- 
sipation range (see next section) might not be sufiiciently 
resolved in MHD simulations [Verscharen et al., 2012]. 

4. Comparison With Spacecraft Data 

Analytically, particle scattering mean free paths in the 
interplanetary system can be obtained by fitting cosmic-ray 
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Figure 2. The electron meam free path as a function of 
the particle rigidity for different dissipation properties. 
Whereas the black diamonds show the same case as in 
Fig. 1, the other results use fed = 3 x 10* (AU)~^ together 
with p = 3 (gray triangles) and fed = 3 x 10^ (AU)^^ 
together with p = 3.5 (white circles). For clarity, the hor- 
izontal positions are slightly shifted against each other. 

observations to diffusion models. Palmer [1982] concluded 
that, for rigidities between 0.5 and 5000 MV, the parallel 
mean free path is 0.08 AU < Ay < 0.3 AU. This regime is 
usually visualized by the Palmer consensus range [Bieber 
et at., 1994]. In Figs. 1 and 2, this box is compared to nu- 
merically calculated parallel mean free paths by using the 
best-fit parameters in the aforementioned Padian code. In 
Fig. 1, the mean free paths are shown for electrons and for 
protons as the particle rigidity varies. 

In Fig. 1, the original measurements arc shown in addi- 
tion to the Palmer consensus range. Such is done because 
the latter was strongly criticized by Reames [1999] because 
it was obtained by excluding impulsive, scatter free events 
and Eastern solar events, thereby leading to values of the 
parallel mean free path smaller than the general case could 
be. Indeed, both the numerical and the experimental results 
shown in Fig. 1 correspond to mean free paths substantially 
larger than those of the Palmer consensus. It is interest- 
ing to notice that the numerical results compare well with 
the experimental data, even beyond the range of validity of 
Palmer consensus. 

In Fig. 2, the electron mean free path is shown for vary- 
ing dissipation range parameters ka and p. Thereby, the 
effects of the dissipation range and the dissipation spectral 
index are illustrated. For small rigidities, the mean free 
path is generally increased if the dissipation wavonumbcr, 
fed, and/or the dissipation spectral range, p, arc increased. 
Furthermore, Fig. 2 shows that the electron mean free path 
can be as large as 2 AU, which is reminiscent of the works 
by Perri and Zirnbardo [2007, 2008], who inferred diverg- 
ing mean free paths for electrons in the solar wind. The 
importance of alternative acceleration mechanisms without 
pitch-angle scattering is emphasized, which however can- 
not be reproduced in the simulations presented here be- 
cause the magnetic turbulence model always leads to pitch- 
angle scattering for all rigidity values considered. Especially 
a superdiffusive mean free path [Tautz and Shalchi, 2010] 
requires either purely two-dimensional turbulence [Shalchi 
et al., 2008] or significantly stronger electric fluctuations 
[e.g., Tautz, 2010b]. Also, electron scatter free events are 
often described in the literature [e. g., Haggerty and Roelof, 
2002; Lin, 2005]. 



Figure 3. The ratio of the perpendicular and paral- 
lel mean free path as a function of the particle rigidity. 
Shown are the results for electrons (black diamonds) and 
for protons (light squares) together with the respective 
error bars. The dashed and dot-dashed lines show the av- 
erage values for electrons and protons, respectively. For 
comparison, the Palmer consensus range [Palmer, 1982; 
Shalchi et al, 2006] is shown. For clarity, the horizontal 
positions are slightly shifted against each other. 

By comparing simulation results and observations, ones 
notes that the simulation data agrees with the bulk of the 
experimental events. In Fig. 1 of Bieber et al. [1994], it can 
be seen that there arc a inimber of electron events around 
0.1 MV with moan free paths of roughly Ay ~ 1 AU, a fea- 
ture that is accurately reproduced by our simulation results. 
Likewise, there arc proton events around lO'^-lO'' MV with 
A|| « 1 AU, which is equally well reproduced by the numeri- 
cal results shown here. Therefore, the excellent agreement in 
comparison to the Palmer consensus range confirms the var 
lidity of the underlying model that has been constructed to 
describe the dynamics of the turbulent magnetic field. Fur- 
thermore, Fig. 3 illustrates the transport anisotropy, which 
is expressed by ratio of the perpendicular and parallel meam 
free paths. The average ratio is X±/X\\ « 0.0379 ± 0.0093 
for electrons and Ax/Ay ~ 0.0597 ± 0.0426 for protons, in 
agreement with the Palmer values [Palmer, 1982; Shalchi 
et al., 2006]. 

It should be noted that, in the numerical simulations, the 
effect of the turbulent electric field has been neglected inten- 
tionally in order to save computation time. Such is justified 
by the fact that, for shear Alfvcn waves, the magnitude of 
the electric field is reduced by a factor of va/c « 8.9 x 10~^ 
as compared to the magnetic field amplitude [Schlickeiser, 
2002]. An additional simulation run confirmed that the in- 
clusion of the turbulent electric field has a negligible effect 
on the spatial diffusion coefficients. 

Furthermore, there has been a strong debate as to the 
form of the spectrum in the dissipation range, which could 
be either a second inertial range [Sahraoui et al., 2009] or a 
true dissipation range [Alexandrova et al., 2009]. However, 
the modeling of the dissipation range in Eqs. (4) and (5) is 
purely heuristic and is therefore independent on the nature 
of that spectral range. Instead, the simulation results are 
sensitive to the features of the "dissipation" range or, more 
precisely, on the parameters chosen for kd and p so that new 
findings might alter the particle behavior. In the literature 
[see Sec. 4.3 of Bieber et al., 1994, and references therein] it 
is argued that the spectral index must ultimately be larger 
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than 3 for theoretical reasons but, at the same time, 3 pro- 
vides a hmit beyond which the mean square magnetic curl 
diverges. Therefore, the dissipation range spectral index was 
chosen as p = 3 [cf. Smith et al., 1990]. 

5. Summary and Conclusion 

A fundamental problem in cosmic-ray physics is the the- 
oretical explanation of observed particle mean free paths in 
the solar system. Palmer [1982] described the conflict be- 
tween observations and standard analytical theories. Bieber 
et al. [1994] revisited the problem by combining the quasi- 
linear transport theory with a more advanced model of inter- 
planetary turbulence. This model takes account of spectral 
anisotropy, dissipation effects, and the turbulence dynam- 
ics. However, the applicability of quasi-linear theory used 
by these authors is questionable [Shalchi, 2009] and is, thus, 
not reliable. 

An alternative approach to obtain scattering parameters 
theoretically is the application of test-particle codes. In such 
simulations there is no need for analytical transport theo- 
ries. However, previous simulations have been performed 
for very simplified models of turbulence (e. g., static turbu- 
lence, no dissipation range in the spectrum). Therefore, no 
matter whether analytical or numerical descriptions of the 
transport had been used, the theory of particle diffusion has 
always been incomplete. 

In the present paper we combined a realistic model for 
solar wind turbulence with the advanced Padian code. Nu- 
merical results for the parallel mean free path Ay were com- 
pared with Palmer's consensus range in Figs. 1 and 2. Fur- 
thermore, Fig. 3 shows the ratio of the perpendicular and 
parallel mean free paths, X±/X\\. According to the compar- 
ison shown here, solar wind observations of energetic parti- 
cles can be reproduced with excellent accuracy. This posi- 
tive result also confirms our present understanding of inter- 
planetary turbulence and, accordingly, the slab/2D model. 

It is important to note that all parameters were chosen 
according to standard interpretations of solar wind observa- 
tions as, e.g., for; (i) the relative turbulence strength; (ii) 
the slab/2D ratio; (iii) the time-dependence of the slab com- 
ponent; (iv) the Alfvcn velocity; (v) the dissipation range 
spectral index and wavenumber. Whereas a full parame- 
ter study was beyond the scope of the present investigation, 
instead we emphasize that the Bieber et al. [1994] observa- 
tions of interplanetary particles were reproduced with excel- 
lent agreement by using test-particle sinmlations. Therefore, 
such theoretical results do no longer depend on uncertain 
analytical transport theories such as quasilinear theory. In 
light of these new results, it can be concluded that the de- 
velopment of analytical theories may still be based on the 
standard parameters and assumptions listed before, if: (i) 
the slab spectrum is chosen according to the new hybrid 
model, and (ii) a dissipation range is included. 
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